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SOME  COMMON -SENSE  OPTIMIZATION  TECHNIQUES  FOR 
NON- DIFFERENTIABLE  FUNCTIONS  OF  SEVERAL  VARIABLES 

ABSTRACT : 

The  problem  of  obtaining  global  optima  of  non-dif ferentiable 
functions  of  several  variables  is  studied.  In  general,  the  func¬ 
tions  are  multimodal  and  continuous  on  a  compact  domain.  Two 
distinct  methods  are  proposed  and  to  some  extent  compared:  The 
method  of  systematic  search  and  the  random  search  technique.  The 
method  of  uniform  saturation  [the  one  variable  version  of  the 
systematic  search  method]  is  based  on  bisecting  the  interval  (in 
the  one-variable  case)  repeatedly.  Without  loss  of  generality, 
we  may  restrict  the  discussion  to  the  closed  unit  interval 
I  =  [0,1].  At  the  first  stage,  n  =  1,  bisect  the  interval  I  using 
the  point  x  =  1/2.  Let  Mj^  =  max  [f  (1/2)  ,f  (1)  )  .  At  the  second 
stage,  n  =  2,  bisect  each  of  the  intervals  [0,1/2]  and  [1/2,1]  us¬ 
ing  the  points  x  =  1/4  and  x  =  3/4  respectively.  Let  M2  = 

[m^,f (1/4) ,f (3/4) ] .  By  the  nth  stage  we  would  have  subdivided 
the  interval  I  into  2n  subintervals,  each  of  length  (1/2) n,  where¬ 
in  the  partition  points  over  and  above  those  previous  stages  are 
i(l/2)n,  i  =  1,3, . . . ,2n-l.  Thus  the  Mn's  are  inductively  given  by 

Mn  =  max [Mn-i , f (i/2n) ;  i  =  1 , 3  , . . . 2n— 1 ] .  It  is  now  clear  that  Mn 
is  monotonic  increasing  sequence,  <_  M2  £  M2  £  .  .  .  If  we  re¬ 
peat  the  procedure  enough  times,  we  would  "saturate"  the  interval 
I  by  evenly  spaced  points  in  such  a  way  that  the  distance  between 


two  neighboring  points  diminishes  geometrically  as  n-increases. 

Thus  we  "zero-in"  on  a  solution  of  the  problem.  This  method  is 
later  modified  to  the  case  of  functions  of  two  or  three  variables. 

The  Random  Search  Technique  used  here  determines  all  the 
optimal  points  of  the  non-dif ferentiable  continuous  functions  with 
many  variables  defined  on  compact  domain.  The  procedure  begins 
with  evaluating  the  given  function  at  pre-determined  number  of 
points  selected  randomly  over  the  closed  bounded  domain.  Suppose 
m  points  are  selected  randomly  over  the  domain  and  the  function 
is  evaluated  at  each  of  the  m  points.  The  minimum  functional 
value  and  the  point  at  which  the  minimum  occurs  (if  the  problem  is 
one  of  minimization)  are  saved.  This  step  is  carried  out  n 
times,  where  n  is  sufficiently  large.  The  resulting  n  points 
will  cluster  around  the  minima.  Suppose  there  are  r  cluster 
points,  then  there  is  a  possibility  that  around  each  cluster  point, 
a  local  minimum  may  exist.  We  develop  a  single  program  to  find 
all  the  cluster  groups  as  well  as  cluster  points  using  a  local 
optimization  routine.  Thus  the  global  minimum  is  obtained  by 
simple  comparison.  The  new  method  developed  here  is  clearly  an 
improvement  with  regards  to  time  and  accuracy  over  the  methods  pro- 
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INTRODUCTION 


There  are  many  optimization  procedures  which  enable  one  to 
determine  the  minimum  of  a  unimodal  function  in  n-space.  If  the 
function  is  differentiable  in  a  compact  domain,  global  minima  may 
be  obtained  through  the  use  of  derivatives.  However,  the  pro¬ 
blem  of  global  optimization  of  multimodal  function  has  received 
comparatively  little  attention,  more  so  when  the  function  in  ques¬ 
tion  is  non-dif ferentiable .  No  efficient  method  has  been  developed 
to  tackle  glrbal  optimization  problems. 

As  a  general  principle,  the  accuracy  with  which  a  procedure 
locates  optima  improves  with  the  number  of  functional  evaluations. 
In  principle,  however,  one  seeks  a  balance  between  a  degree  of 
certainty  and  the  cost  of  implementation.  A  procedure  which  lo¬ 
cates  optima  with  great  precision  and  certainty  would  be  practical¬ 
ly  worthless  if  it  requires  economically  unfeasible  number  of  cal¬ 
culations. 

There  are  several  methods  presently  utilized  to  seek  global 
optimua;  among  them  are  those  suggested  by  Brooks  [1] ,  Becker  and 
Lago  [2]  and  Price's  CRS  method  [3].  The  Simple  Random  Method  ac¬ 
cepts  the  optimum  function  value  as  global  optimum  after  making  a 
specified  number  of  trials  randomly  selected  from  the  domain.  The 
stratified  Random  Search  method  divides  the  domain  into  a  number 
of  subdomains  of  equal  size  and  selects,  at  random,  a  trial  point 
from  each  subdomain  and  each  time  keeps  the  optimal  function  value. 
The  procedure  is  repeated  a  good  many  times.  Some  improvement  on 
the  simple  random  search  is  provided  by  Becker  and  Lago.  Their 


procedure  begins  with  a  Simple  Random  Search  over  the  domain,  in¬ 
stead  of  retaining  the  single  point  with  the  optimal  function 
value,  Becker  and  Lago  retain  a  predetermined  number  of  points  with 
optimal  function  values  in  each  trial.  If  the  number  of  trials  is 
sufficiently  high,  the  retained  points  tend  to  cluster  around  some 
optima.  Then  a  mode  seeking  algorithm  is  used  to  group  the  points 
into  discrete  clusters  and  to  define  the  boundaries  of  the  sub- 
regions  each  embracing  a  cluster.  The  clusters  are  graded,  by 
searching  in  each  for  the  retained  points  with  the  lowest  function 
value  and  then  rated  according  to  the  relative  values  of  the  clus¬ 
ter  minima.  The  entire  procedure  is  then  repeated  using  as  the 
initial  search  region  that  subdomain,  defined  by  the  mode  seeking 
algorithm  around  the  'best'  cluster.  The  user  may  choose  to  ex¬ 
amine  also  the  second  best  cluster,  or  indeed  all  clusters,  accord¬ 
ing  to  the  extent  of  his  doubt  as  to  whether  or  not  the  global 
minimum  will  be  found  in  the  subdomain  defined  by  the  best  cluster. 

The  controlled  Random  Search  (CRS)  suggested  by  Price  is 
similar  to  Becker  and  Lago,  but  CRS  combines  the  random  search  and 
mode-seeking  algorithm  into  a  single  continuous  process.  But  the 
problems  of  inefficiency  and  economic  consideration  still  remain. 

METHOD  OF  SOLUTIONS 

This  paper  deals  with  two  methods:  (I)  Systematic  Search  (The 
Method  of  Uniform  Saturation),  (II)  Random  Search.  ^In  both  cases 
it  is  assumed  that  the  functions  are  defined  and  continuous  on  a 
compact  domain.  They  are  also  assumed  to  be  multimodal  functions. 

In  general  the  systematic  search  does  not  provide  all  the  optimal 
points,  the  primary  emphasis  here  being  location  of  a  global  optimum. 


Despite  several  restrictions  and  difficulties,  the  Random  Search 
method  attempts  to  obtain  all  the  optima,  one  optimum  point  in 
each  mode . 

I.  SYSTEMATIC  SEARCH  (The  Case  of  Two  or  Three  Variables) 

Suppose  f(x,y)  is  continuous  on  the  closed  unit  square 
S  =  {(x,y):0  _<  x,y  <  1}.  Then  certainly  a  subdivision  of  the  inter¬ 
val  [0,1]  into,  say,  50  equal  subintervals  would  have  to  be  con¬ 
sidered  as  a  reasonable  partition.  That  is  to  say,  50  is  a  reason¬ 
ably  small  number.  Yet  even  with  50  partition  points  on  each  of 
the  x  and  y  axes,  we  are  faced  with  50  x  50  =  2500  partition 
points  of  the  unit  square  S.  For  the  case  of  a  function  f(x)  of 
one  variable,  we  certainly  would  want  to  partition  the  interval 
[0,1]  into  MORE  than  50  subintervals  to  get  a  reasonable  assurance 
that  an  optimum  has  been  included.  Hence  we  cannot  be  confident 
that  the  global  optimum  will  be  among  the  values  at  the  2500  parti¬ 
tion  points  of  the  square  S. 

The  case  of  a  function  of  three  variables  is  much  worse.  Here 
a  subdivision  of  each  co-ordinate  AXIS  into  50  partition  points  re¬ 
sults  in  50  x  50  x  50  =  125,000  points  of  the  cube.  This  is  just 
to  get  a  crude  starting  point.  Hence  we  see  that  the  number  of 
evaluations  becomes  prohibitive  very  rapidly  and  so,  to  have  any 
hope  whatsoever  of  handling  the  multivariable  case,  we  would  have 
to  abandon  the  purely  exhaustive  scheme  (Method  of  Uniform  Satura¬ 
tion)  used  in  the  uni variable  case. 

The  proposed  method  is  based  on  two  steps.  The  first  step 
involves  consideration  of  an  initial  grid  on  the  domain.  An  initial 


< 

I 


point  is  then  obtained  based  on  the  grid.  The  second  step  starts 
with  the  initial  point  and  proceeds  by  the  method  of  'crossings'. 

Any  direct  search  procedure  such  as  the  one  presently  given 
would  require  a  large  number  of  evaluations.  For  a  function  f(x,y) 
of  two  variables  on  a  rectangle,  we  consider  100  partition  points 
on  each  of  the  x  and  y  axes  to  be  reasonable.  This  gives  rise 
to  100  x  100  =  10,000  partition  points.  We  realize  that  10,000 
evaluations  might  not  be  cost-effective  and  that  other  more  effi¬ 
cient  methods  might  be  employable.  The  fact  remains  that  this 
procedure  is  direct,  simple  to  execute  and  self-contained  (not 
based  on  other  search  procedures  already  in  existence) . 

Several  theorems  pertaining  to  functions  of  two  variables  are 
proved  and  some  twenty  one  illustrative,  computational  examples  are 
provided.  These  examples  comprise  Tables  1,  2  and  3.  The  computer 
programs  are  given  in  Appendix  A. 
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The  One-Variable  Case:  (The  Method  of  Uniform  Saturation) 

Consider  the  non-linear  programming  (NLP)  problem:  MAXIMIZE 

f  (x)  :  a  <_  x  <_  b,  where  f  :  I  +  R  is  a  continuous  real-valued 
function  defined  on  the  closed  interval  I  =  [a,b].  Without  loss  of 
generality,  we  may  restrict  the  discussion  to  the  closed  unit  in¬ 
terval  I  =  [0,1].  At  the  first  stage,  n  =  1,  bisect  the  interval 
I  using  the  point  x  =  1/2.  Let  My  =  max  {f  (1/2)  ,f (1) } .  At  the 
second  stage,  n  =  2,  bisect  each  of  the  intervals  [0,1/2]  and 
[1/2,1]  using  the  points  x  =  1/4  and  x  =  3/4  respectively.  Let 
M2  =  max  My ,f (1/4) ,f (3/4)  .  At  the  third  stage,  n  =  3,  bisect 

each  of  the  intervals  [0 , 1/4 ] , [1/4 , 1/2 ] , [1/2 , 3/4 ] ,  and  [3/4,1]  us¬ 
ing  the  points  x  =  1/8,  x  =  3/8,  x  =  5/8,  and  x  =  7/8  respectively. 
Set  My  =  max{M2 , f (i/8)  :  i  =  1,3, 5, 7}. 

By  the  n'th  stage  we  would  have  subdivided  the  interval  I  in¬ 
to  2n  subintervals,  each  of  length  (1/2) n,  wherein  the  new  parti¬ 
tion  points  over  and  above  those  of  the  previous  stages  are 
i(l/2)n  :  i  =  1,3,. ..,2n  -  1.  Thus  the  Mn's  are  inductively  given 
by  Mjj  =  max  {M^y ,  f  (i/2n)  :  i  =  l,3,...,2n  -  1}.  It  is  now  clear 

that  Mn  is  monotone  increasing,  viz.  My  <.  M2  <_  My  ...  If  we  repeat 
the  procedure  enough  times,  we  would  "saturate"  the  interval  I  by 
evenly  spaced  points  in  such  a  way  that  the  distance  between  two 
neighboring  points  diminishes  geometrically  as  n  increases.  Thus 
we  "zero  in"  on  a  solution  of  the  problem.  That  is,  if  xQ  SOLVES 

the  problem,  then  there  is  a  bisecting  point  x^  WITHIN  ANY  PRE¬ 
SCRIBED  DISTANCE  from  xQ.  Thus  if  e  >  0  is  preassigned,  we  are  as¬ 
sured  of  the  existence  of  an  x^  for  which  | x^  -  xQ |  <  e  whenever 


n  j  -ii '  iiiifnr iij-  iiiiriiir  "jar 


n  is  such  that  2n  >  1/e.  Since  the  function  f (x)  is  continuous, 
we  know  that  ffx^)  will  be  close  to  f (xQ)  whenever  x^  is  "suffi¬ 
ciently"  close  to  xQ. 

The  Two-Variable  Case 

We  next  consider  a  real-valued  function  f(x,y)  which  is  con¬ 
tinuous  on  the  closed  unit  square  S  =  {(x,y):  0  <_  x,y  £  1}.  The 
non-linear  programming  (NLP)  problem  is:  MAXIMIZE  f(x,y)  :  (x,y)  e  S. 

Theorem  1  Let  f(x,y)  be  a  real-valued  function  which  is  continu¬ 
ous  on  a  compact  domain  D.  Then 

MAXIMUM  f(x,y)  =  MAX  {MAX  f(x,y)},  where  for 
(x,y)eD  xeDy  yeDx 

each  fixed  x,  Dx  =  {y  :  (x,y)e  D}  and  for  each  fixed 
y,  Dy  =  x: { (x,y)  e  D). 

Proof  Clearly  MAXIMUM  f(x,y)  >_  MAX  {MAX  f(x,y)}.  Suppose 

(x,y)eD  xeDy  yeDx 

that  the  inequality  is  strict: 

MAXIMUM  f (x,y)  >  MAX  {MAX  f(x,y)}.  Say 
(x,y)eD  xeDy  yeDx 

MAXIMUM  f (x,y)  =  f ( Xq ,  yQ) 

then  f(xQ,  yQ)  >  MAX  {MAX  f(x,y)} 

xeDy  yeDx 

>  MAX  f ( x  ,  y) 

y£Dx 

>  f  (Xq,  yQ)  ,  a 

contradiction . 

As  a  corollary ,  we  have: 

If  f(x,y)  is  continuous  on  the  closed  unit 


square  S ,  then 


MAXIMUM  f (x,y) 
(x,y) eS 


=  MAX  {  MAX  f (x,y) }  = 
0<x_<  1  O^yjcl 

MAX  {MAX  f (x,y) } . 

0<y£l  0<x_<l 


We  point  out  that  the  assumption  of  continuity  cannot  de 
weakened  to  separate  continuity  as  the  following  example  shows. 


Example  1 

f (x,y) 


fxy/(x4  +  y4)  if  (x,y)  e  S  -  (0,0) 
0  if  (x,y)  =  (0,0) 


Theorem  2  Let  f(x,y)  be  continuous  on  the  unit  square  S.  For 

each  a  e  [0,1],  define  h(a)  =  MAXIMUM  f(a,y).  Then  the 
function  h:  [0,1]  +  R  is  continuous. 


Proof 

Let  a  e  [0,1]  and  let  e  >  0  be  given.  Uniform  con¬ 
tinuity  of  f(x,y)  implies  the  existence  of  6  =  6(e)  >  0 
such  that  |f(a,y)  -  f (x,y) |<  e  whenever  |x  -  a|  <  6. 

Take  x  such  that  |x  -  a|  <6  and  let  the  maximum  of 

f(a,y)  over  y  occur  at  y  and  let  the  maximum  of  f(x,y) 

over  y  occur  at  y.  That  is,  h(a)  =  MAXIMUM  f(a,y)  = 

0<y<L 

f(a,y)  and  h(x)  =  MAXIMUM  f(x,y)  =  f(x,y). 

0<y<l 

Then  jf(a,y)  -  f(x,y)|  <  e  and  |f(a,y)  -  f(x,y)|  <  e 
f(a,y)  <  f(x,y)  +  e  <_  f(x,y)  +  e  and  f(a,y)  -  f(x,y)>-e 

f(a,y)  -  f(x,y)  <  c  and  f(a,y)  -  f(x,y)>-e. 

Thus  |f(a,y)  -  f(x,y)|  <  e  whenever  |x  -  a|  <  6  or 
|h(a)  -  h(x)  j  <  e  whenever  |x  -  a|  <  6.  This  shows 
h:  [0,1]  ■»  R  is  uniformly  continuous  on  [0,1]. 


The  example  cited  earlier  shows  that  the  assumption  of  con¬ 
tinuity  on  f(x,y)  in  Theorem  2  cannot  be  relaxed  to  separate  con¬ 


tinuity. 


Example  2 


f(x,y)  = 


fxy/(x4  +  y4)  if  (x,y)  e  S  -  {(0,0)} 
0  if  (x,y)  =  (0,0). 


It  is  easily  verified  that  here  the  function  h(a)  =  MAXIMUM  f(a,y) 

0<y<l 

is  given  by: 


h  (a) 


( 


3 

4  <TT?> 


0 


if  0  <  a  <  1 


if  0  =  a.  That  is,  h(a)  is  not 


continuous  at  a  =  0. 

The  next  theorem  appeared  as  Problem  E  2854  and  its  solution 
in  the  April  1982  issue  of  the  American  Mathematical  Monthly. 


Theorem  3  Let  f(x,y)  be  a  real-valued  continuous  function  on  the 

unit  square  S  =  {  (x,y)  :  0  x,y  <_  1}.  Additionally, 

suppose  that  for  each  a  e  [0,1],  the  maximum  of  f(a,y) 

over  y  occurs  at  ONLY  ONE  value  of  y,  say  MAXIMUM  f(a,y) 

°<y<_l 

f(a,y*(a)).  Then  the  assignment  a  b  y* (a)  defines  a 
continuous  function  y*  :  [0,1]  -+■  [0,1]. 

Proof  The  proof  is  by  contradiction.  Suppose  y*  is  not  con¬ 

tinuous  at  some  a  e  [0,1].  Let  {an}  be  a  sequence  in 

[0,1]  convergent  to  a  such  that  bn  =  y*(an)  fails  to 
converge  to  y  (a) .  Since  { bn }  is  a  sequence  in  a  com¬ 
pact  space,  we  may  assume,  without  loss  of  generality, 


that  {bn}  converges  to  some  number  b  e  [0,1]  (otherwise  we  select 

a  convergent  subsequence).  Let  f(a,y*(a))  -  f(a,b)  =  e.  Since 

f(a,y*(a))  =  MAXIMUM  f(a,y)  we  see  that  f(a,y*(a))  _>  f(a,b);  i.e.  , 
0<y<l 

e  ^  0.  The  uniqueness  of  the  maximum  implies  that  e  >  0.  For  if 
e  =  0  then  f(a,y*(a))  =  f(a,b)  or  BOTH  y*(a)  AND  b  maximize  f(a,y) 
and  y*(a)  =  b,  violating  the  assumed  uniqueness  of  the  maximum. 

We  have: 

jf (a,y*(a))-f (a,b) |  <  | f (a ,y* (a) ) -f (an ,bn) |  +  j f (an ,bn) -f (a,b) | 

=  |h(a)  -  h(an) |  +  | f (an,bn) -f (a,b) | 

where  h  is  as  defined  in  Theorem  2. 

Theorem  2  together  with  the  fact  that  an  a  implies  that 

|h(a)  -  h(an) |  <  l/2e  whenever  n  is  sufficiently  large.  Also,  con¬ 
tinuity  of  f(x,y)  together  with  the  convergences  aR  a  and  bn  -*•  b 
implies  |f(an,bn)  -  f(a,b) |  <  l/2e  whenever  n  is  sufficiently  large 

Thus  taking  n  so  large  that  BOTH  l/2e-inequalities  hold  simultane¬ 
ously  we  obtain  the  following  contradiction: 

f (a,y* (a) )  -  f (a,b) |  <  l/2e  +  l/2e 

e  <  e . 

We  acknowledge  our  gratitude  to  Dr.  Charles  Giel  (formerly  of  A&T 
State  University)  for  the  proof  of  Theorem  3  above. 

In  a  private  communi fi cation ,  Professor  R.  A.  Struble  of  North 
Carolina  State  University,  gave  the  following  solution  to  Problem 
E  2854  and  hence  an  independent  proof  of  Theorem  3. 

Alternate  Proof  of  Theorem  3  (Direct  Proof) 

Let  a  [0,1]  be  given  and  let  {an}  be  a  sequence  in  [0,1]  such 

an  *  a*  We  sllow  Y*<an)  Y*(a).  The  sequence  {y*(an)}  is  in  the 
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compact  space  [0,1]  and  hence  we  may  assume  that  (y*(an)}  is  con¬ 
vergent  to  some  number  be  [0,1]  (otherwise  we  select  a  convergent 
subsequence).  Continuity  of  f(x,y)  implies  that  f(an,y*(an))  -*-f(a,b) 
and  f(an,y*(a))  f(a,y*(a)).  From  the  definition  by  y* ,  it  fol¬ 

lows  f (a  ,y*(an))  >  f(an,y*(a)).  Thus 

lim  f(an,y*(an))  >_  lim  f(an,y*(a))  or  f(a,b)  f(a,y*(a)).  The 

last  inequality  says  b  maximizes  f(a,y)  over  y  so  that  uniqueness 
of  the  maximum  now  implies  b  =  y*(a);  i.e.  ,  y*(an)  -*■  y*(a). 

The  proof  of  Theorem  3  published  in  the  American  Mathematical 
Monthly  is  shorter  than  either  of  the  proofs  given  here;  however 
the  published  proof  relies  on  a  compact  graph  theorem  and,  in  our 
opinion  is  less  instructive.  Problem  E  2854  asks  if  Theorem  3 
may  be  generalized  as  follows.  Suppose  the  requirement  of  the 
uniqueness  of  the  maximum  is  no  longer  imposed  and  the  function 
y*  :  [0,1]  -*•  [0,1]  is  modified  so  that  y*(a)  =  MIN  (y:y  maximizes 
f(a,y)}.  Does  the  assignment  a  •*  y*  (a)  define  a  continuous  func¬ 
tion  y*  ;  [0,1]  -*■  [0,1]?  The  answer  is  NO!  The  following  counter¬ 
example  is  given  in  the  American  Mathematical  Monthly. 

Example  3 

f(x,y)  =  (x-1/2) (y-1/2 ) .  Here  y  (a)= 

Professor  J.  G.  Mauldron  of  Amherst  College  points  out  that 
the  function  of  Example  3  is  unsatisfactory  because  it  fails  to 
satisfy  the  uniqueness  property  miserably  at  a  =  1/2  in  the  sense 
that  the  set  {y:y  maximizes  f(l/2,y)}  =  [0,1]  and  offers  the  follow¬ 
ing  example  instead. 

Example  4 

f(x,y)  =  (x-y)2.  Here  y*  (a)  = 


1  :  a  <  1/2 
0  :  a  >  1/2. 


0  :  a  <  1/2 
1  :  a  >  1/2 
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For  the  function  f(x,y)  of  Example  4,  the  departure  from  the 
uniqueness  condition  is  MINIMAL  in  the  sense  that  the  {y:y  maximizes 
f(a*y))  is  a  singleton  for  a  4  1/2#  while  the  set  {y:y  maximizes 
f (1/2 ,y) }  =  {0,1}. 

Professor  Mauldron  offers  the  following  example  to  illustrate 
that  the  continuity  requirement  on  f(x,y)  in  Theorem  3  cannot  be 
relaxed  to  separate  continuity. 


Example  5 


/  7  i 

Cx,y)  =  i  2 

|^8y(x-y)/xz  i 


f  x  =  0 


f  x  4  o 


The  function  f(x,y)  satisfies  the  uniqueness  condition  but  is  only 
separately  continuous.  The  induced  function  y*(a)  is  discontinuous 


at  a  =  0: 


*  f  1  if  a  =  0 

'  (a)  =  4 

L  l/2a  if  0  <  a  _<  1 


Looking  at  Examples  3  and  4#  one  may  be  tempted  to  conjecture  that 
y*  :  [0,1]  -*■  [0,1]  enjoys  the  property  of  one-sided  continuity. 
Professor  Richard  Tucker  of  A&T  State  University  gives  the  follow¬ 
ing  counter-example. 


Example  6 


F(x,y)  = 


f(x,y)  :  0  <  x  <_  1/2,  0  <_  y  <  1 
f  (1-x  ,y )  :  1/2  <x<l,  0_<y<l 


where  f(x,y)  is  as  in  Example  4  or  f(x,y)  =  jx  -  y|. 


Here  y  (a)  = 


0  <  a  <  1/2 
a  =  1/2 


1/2  <  a  <  1. 


COMPUTATIONAL  EXAMPLES 


Aaron  Chew  wrote  the  BASIC  programs  for  use  on  Texas  Instru¬ 
ments  99/4A  personal  computer  with  Extended  Basic  module  and 
Peripheral  Expansion  System.  We  express  our  deep  appreciation  to 
Aaron  for  his  programming  assistance. 

The  TWO-VARIABLE  PROGRAM  is  based  on  the  following  procedure. 
Let  f(x,y)  be  defined  on  the  closed  rectangle  R  =  {  (x,y)  :a  x  b; 
c  <  y  <  d  .  First  use  an  Initial  Grid  on  the  rectangle  obtained  by 
putting  evenly  spaced  points  on  the  sides  of  the  rectangle  lying  on 
the  co-ordinate  axes : 

a  =  XQ  <  xx  <  x2  <  . . .  <  xM  =  b;  c  =  yc  <  yx  <  y2  <  . . .  <  yM  =  d 

where  x.  =  a  +  i  (b-a)  and  y .  =  c  +  j  The  procedure  first 

1  M  3  M 

produces  an  initial  approximation  (x,y)  based  on  the  points 
(Xi,yj)  of  the  initial  grid.  The  Main  Program  then  uses  (x,y)  as 

STARTING  POINT  and  procedes  as  follows.  Fix  x  =  x  and  minimize 
f(x,y)  over  y  e  [c,d]  using  evenly  spaced  partition  of  the  type 

N 

used  in  the  one-variable  case;  namely,  evenly  spaced  points  (1/2) 

apart.  Say  MIN  f(x,y)  occurs  at  y  =  y  .  Next  minimize  f(x,y  ) 

Y  11 

over  x  e  [a,b] ,  again  using  points  that  are  (1/2)N  apart.  Say 

MIN  f(x,y_)  occurs  at  x  *  x,.  Refer  to  (x_,y»)  as  the  second 

x  2  2  2  2 

CROSSING.  Repeat  as  often  as  desired. 
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II.  RANDOM  SEARCH 

The  domain  of  the  function  is  closed  and  bounded  and  it  will 
always  be  possible  to  select  the  initial  starting  points  at  the 
boundary.  All  the  examples  discussed  here  are  of  functions  whose 
domains  are  of  the  shape  of  hypercubes ,  a^  <  <  b^.  Therefore, 

starting  points  may  be  taken  as  a^ ,  i  =  1,2,3,...  The  next  point 
may  be  taken  as  a^  +  e,  where  e  =  (b^  -  a  )/N,  if  one  decides  to 

use  N  points  to  obtain  the  first  minimum.  It  is  not  really  impor¬ 
tant  which  formula  is  used  to  generate  points  over  the  domain,  as 
long  as  those  domains  are  searched  repeatedly  without  duplication. 

We  evaluate  at  the  first  N  points  just  generated  and  store  the 
minimum  and  the  coordinates  of  the  minimizing  point.  We  repeat  the 
procedure  M  times.  Therefore,  in  all  M  minimum  values  are  saved 
together  with  the  coordinates  of  the  minimizing  points.  All  the 
generated  points  have  to  be  tested  whether  they  belong  to  the  do¬ 
main  before  they  can  be  used.  The  essential  features  of  the  al¬ 
gorithm  are  indicated  in  the  f low-diagram  (Figure  1) . 

The  M  stored  points  should  cluster  around  the  minima.  An  il¬ 
lustration  of  this  concept  is  shown  in  Figure  2.  The  main  task  of 
this  procedure  is  to  locate  all  the  cluster  groups.  We  have  achiev¬ 
ed  only  partial  success  in  reaching  this  objective  because  of  a  pro¬ 
blem  described  below: 

If  some  of  minima  lie  very  near  to  each  other,  this  procedure 
cannot  separate  the  clusters,  because  the  radius  of  the  hypersphere 
which  embrace  these  cluster  points  should  be  very  small  and  there¬ 
fore  many  points  still  remain  outside  of  any  hypersphere.  These 
points  which  are  outside  give  false  cluster  groups  and  thereby  in¬ 
crease  the  function  evaluations  later  tremendously.  Let  us  take  the 


EVALUATE  fC*i)  at 
Ac;  G ENER/ATE  U-  M0R£ 
TOUTS 


FIGURE  1 


following  function  as  example:  f(x^,x2)  =  (|x1 1— 1) 2  +  (|x2|-2)2 
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-  4  <  *1»x2  <_  4 . 

There  are  four  minima  with  function  value  f  =  0  and  coordinates 
(-1,-2) , (-1,2) , (1,-2) , (1 ,2)  .  All  these  clusters  lie  quite  far  a- 
part  and  our  procedure  can  obtain  all  of  them  very  quickly.  How¬ 
ever,  if  f(x1,x2)  =  ( |  x-^  |  -  0.1) 2  +  (|x2|  -  0 . 5 ) 2  this  procedure 

will  lead  to  one  minimum  point  only  since  all  four  points  (0.1,0. 5), 
(0. 1,-0. 5)  ,  (-0.1, 0.5)  and  (-0.1, -0.5)  are  lying  on  a  very  small 
rectangle . 

After  separating  the  cluster  the  next  biggest  task  is  to  find 
the  actual  minimum  in  each  cluster  group.  Any  local  optimizing 
method  may  be  used.  However,  Nelder  and  Mead  Simplex  Search  method 
is  the  most  efficient  one  for  non-dif ferentiable  functions.  We  have 
used  Nelder  and  Mead  Simplex  Search  method  [4]  in  our  program.  The 
Nelder  and  Mead  Simplex  Search  requires  m  +  1  points  for  m-dimension- 
al  space  and  they  may  not  lie  on  the  same  hyperplane.  Therefore, 
each  cluster  group,  or  the  hyperspheres  which  embrace  the  cluster 
groups  must  include  at  least  m  +  1  points  to  start  the  initial 
simplex.  So,  not  only  is  the  counting  of  points  necessary  in  each 
cluster  but  also  sometimes  the  points  must  be  regenerated  if  the 
points  fall  short. 

The  checking  of  collinearity  is  another  important  task  in 
Simplex  Search  method.  If  the  simplex  repeats  itself  for  a  specific 
number  of  times,  this  has  to  be  modified  to  prevent  from  collaps¬ 
ing  the  simplex.  One  way  to  solve  this  problem  is  to  replace  a 
point  from  the  collapsing  simplex  by  a  point  which  lies  on  the  or¬ 
thogonal  direction  to  the  hyperplane. 


‘\D  -■ 

r  V 
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The  Choice  of  the  Number  of  Retained  Points; 

The  number  of  retained  points  may  depend  on  the  size  of  the 
domain  and  as  well  as  the  number  of  variables.  Suppose  we  start 
with  fifty  points;  fifty  functional  evaluations  are  performed  and 
one  point  with  lowest  functional  value  is  retained.  If  one  wants 
to  retain  50  points,  2500  function  evaluations  are  required.  There¬ 
fore,  the  number  of  function  evaluations  is  very  high  where  as  stor¬ 
age  requirement  is  comparatively  less. 

Constraints ; 

All  the  global  optimization  problems  may  be  regarded  as  con¬ 
strained  in  the  sense  that  the  search  is  confined  within  the  initi¬ 
ally  prescribed  domain.  If  any  point  falls  out  of  this  domain, 
that  point  has  to  be  discarded.  When  additional  constraints  are 
imposed,  then,  depending  on  the  number  and  complexity  of  these  con¬ 
straints,  a  sufficiently  large  number  of  points  has  to  be  selected 
to  insure  that  a  reasonable  proportion  of  points  from  the  totality 
of  trial  points  be  included. 

The  program  is  written  in  FORTRAN  IV  and  several  examples  are 
discussed.  Since  we  are  using  Simplex  Search  Method,  the  number  of 
dimensions  must  be  more  than  one.  The  program  is  attached  in  Ap¬ 
pendix  B. 

Example  1 

The  function  to  be  minimized  is 
f(x,y,z)  =  (x  -  y  t  z)2  +  (-x  +y  +z)2 

+  (x  +  y  -  z)2, 

-  1  <_  x,  y,  z  <_  1 


It  is  easy  to  show  that  f  is  a  strictly  convex  quadratic  function 
with  an  unique  minimum  at  (0,0,0)  and  f  =  0. 


After  2732  iterations,  we  have 
f  =  0 
x  =  0 
y  =  0 
2  =  0 

Actual  values:  f=0,x=0,y=0,  z=0 

The  method  of  systematic  search  takes  12096  iterations  to  arrive 
at  this  result. 

In  this  connection  it  must  be  pointed  out  that  in  using  the 
systematic  search  method,  we  have  tried  to  adhere  to  standardized 
values  for  the  number  of  initial  grid  points  and  the  number  of 
crossings.  Since  the  function  is  NON-NEGATIVE  and  the  actual  opti¬ 
mal  point  is  (0,0,0)  ,  the  method  of  bisection  would  yield  the 
answer  on  the  very  first  bisection  (27  evaluations  at  most!).  Hence 
the  computer  operator  would  STOP  the  computer  after  ONLY  27  evalu¬ 
ations  because  he  sees  that  f  already  attains  0  [and  can  never  be 
improved]  after  27  evaluations. 

Example  2 

This  example  is  used  to  compare  the  result  obtained  by  the 
method  systematic  search  (discussed  in  this  paper,  Example  19, 

Table  3)  and  the  actual  values.  The  function  is 
f(x,y,z)  =  I x— 1 |  +  | y-1 . 5 |  +  |6z-l|. 

0  <  x,  y  <  3,  0  <  z  <  1.5. 

Actual  solution: 

Min(f(x,y,z)  =  0 
x  =  1,  y  =  1.5,  z  =  0.1666... 


By  the  Systematic  Search  Method: 
Min(f(x,y,z)  =  0.00390625 
x  =  1,  y  =  1.5,  z  =  0.16605625 
Number  of  evaluations:  18752 


By  the  Random  Search  Method: 

Min(f(x,y,z)  =  0.000001326 
x  =  1,  y  =  1.5,  z  =  0.166667 
Number  of  evaluations:  3008 
Example  3 

As  another  example,  let  us  take  the  following  function  which 
was  chosen  by  both  Becker  and  Lago  and  Price's  CRS  algorithm  (with 
additional  constraint) : 

2  2 

f(x1,x2,x3)  =  9-8x^-6x2“4x2+2x^+2x2 
2 

+X3+2x1x2+2x1x3 

0  1  xlx2l3'  0  1  x3£l.5. 

The  actual  solution  is 

f  =  0,  x^=  1,  x2=  1,  x3=  1.  The  Random  Search  method  achieves 

this  solution  in  2686  evaluations  where  f  =  -0.1192xl(,’  ^ 
xj  =  1,  x2  =  0.9999,  x3  =  1 

The  method  of  systematic  search  takes  14144  evaluations. 

Example  4 

As  a  final  example,  we  like  to  consider  the  following  func¬ 
tion  to  obtain  all  the  four  minima.  Becker  and  Lago  and  Price  also 
discussed  a  similar  function.  Their  function  was 
f(Xl,x2)  =  (IxJ  -  5) 2  +  <|x2|-  5) 2 . 

Price  obtained  all  the  four  minima  around  0(10“6)  after  5000 
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evaluations  but  not  obtained  the  coordinates.  We  take 
f<xi,x2,x3)  =  { IxjJ  -  5) 2  +  (|x2j  -  5) 2 

+  (x3  -  l)2 
-10  <_  x^#x2,x3  <_  10. 

All  the  four  minima  are  obtained  after  4010  evaluations: 


Function  Value 

Coordinates 

0. 8298xl0“10 

X1 

x2 

x3 

0 . 244xl0“9 

5.0 

5.0 

1.0 

0.1591xl0~9 

-5.0 

-5.0 

1.0 

0.1699xl0~9 

-5.0 

5.0 

1.0 

5.0 

-5.0 

1.0 

•Actual  minimum  is  of  course  0  at  all  these  four  points. 

The  computer  printout  of  the  unified  program  is  enclosed  in 
the  Appendix  B. 

Conclusion : 

The  Random  Search  Method  described  in  this  paper  is  not  really 
a  Random  Search.  Besides  the  initial  point  -  generation  technique, 
everything  later  becomes  more  systematic  than  random.  The  method 
seems  to  be  very  efficient  for  problems  wherever  the  Nelder  and 
Mead  Simplex  Search  method  applies.  It  suffers  a  serious  setback 
if  some  of  the  minima  are  very  'close'  to  each  other.  How  close  is 
very  'close'?  This  is  an  open  question.  One  may  use  different  lo¬ 
cal  optimization  techniques  to  avoid  this  situation.  The  problem 
of  collapsing  simplex  may  be  handled  as  suggested  in  this  paper. 

Examining  the  program,  one  discovers  that  the  storage  require¬ 
ment  is  not  as  great  as  it  first  appears.  All  the  initial  points 
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generated  need  not  be  saved.  We  need  only  to  save  the  number  of 
retained  points  which  actually  form  clusters. 

The  work  on  the  Method  of  Systematic  Search  has  generated  some 
mathematical  theory  and,  it  appears  that  more  theoretical  develop¬ 
ments  may  be  possible.  Some  of  the  advantages  of  the  systematic 
search  method  are: 

(1)  it  is  applicable  to  functions  of  a  single  variable 

(2)  it  is  direct 

(3)  it  is  easy  to  execute  (the  problems  of  simplicial  collapse, 

etc.  do  not  appear) 

(4)  it  is  independent  of  other  search  procedures  already  in  exis¬ 

tence 

(5)  it  goes  after  the  global  optimum  without  first  calculating 

local  optima 

(6)  it  is  not  sensitive  to  'nearness'  of  the  local  optima  to  each 

other 

The  method  suffers  from  the  standpoint  of  being  computational¬ 
ly  uneconomical  in  that  the  number  of  evaluations  increases  geome¬ 
trically  with  an  increase  in  the  number  of  variables.  Also  the 
method  of  systematic  search  does  not,  in  general,  obtain  all  the 
local  minima.  This,  in  turn,  may  lead  to  some  doubt  as  to  where 
the  actual  global  minimum  occurs. 

This  is  a  serious  problem  attributed  to  all  procedures  which 
find  global  minimum  without  calculating  derivatives  such  as  the 
method  of  Becker  and  Lago  and  Price's  Controlled  Random  Search  Pro¬ 
cedure.  However,  the  method  of  Random  Search  appears  to  overcome 
this  problem. 

Summary  of  Most  Important  Results: 

(a)  The  following  three  theorems  have  been  established: 

(1)  Theorem  1;  Let  f(x,y)  be  a  real  valued  function  which  is 


Continuous  on  a  compact  domain  D.  Then 

Max  f(x,y)  =  Max  \  Max  f(x,y)  }  where  for  each  fixed  x, 
(x,y)eD  xeDy  £  yeDx  ) 

Dx  =  {y:(x,y)eD}  and  for  each  fixed  y,  Dy  =  {x:(x,y)eD}. 

(2)  Theorem  2 :  Let  f(x,y)  be  continuous  on  the  unit  square  S.  For 

each  a  e  [0,1],  define  h(a)  =  Max  f(a,y).  Then  the  function  h: 

0<y<l 

[0,1]  ■*  R  is  continuous. 

(3)  Theorem  3*:  Let  f(x,y)  be  a  real  valued  continuous  function  on 
the  unit  square  S  =  {  (x,y)  :0  <_  x,y  £  1}.  Additionally  suppose  that 
for  each  a  e  [0,1],  the  maximum  of  f(a,y)  over  y  occurs  at  only 
one  value  of  y,  say  Max  f(a,y)  =  f(a,y*(a)).  Then  the  assignment 

a  | - >y*(a)  defines  a  continuous  function  y*:[0,l]  -►  [0,1]. 

(b)  A  complete  program  to  find  the  various  cluster  groups  of  a 
multimodal  non-dif ferentiable  continuous  function  defined  on  a  com¬ 
pact  domain  and  to  pinpoint  the  minimum  value  of  the  function  at 
each  cluster  group  using  local  optimizing  technique  is  written. 
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function  (f)  DOMAIN  COMPUTED  SOLUTION _ ACTUAL  SOLUTION 

~~~  —  .3  <  x  <  0  r  -  58.392301*13  f  -  58.392305 
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PROGRAM  ONE  (VARIABLE)  (See  Example  #1) 


1  CONS  =  l.E  20 
10  INPUT  "N"  :  N 

20  INPUT  "START  &  END  POINT"  :  B,E 
30  FOR  X  =  B  TO  E  STEP  .  5aN 

40  A  =  100*  (X  -X  A  2)  A  2  +  [6.4MX-.5)  A  2  -  X  -  .6]/\2  ::  GOSUB  100 
50  NEXT  X 
55  GOTO  110 

100  IF  A  <  CONS  THEN  CONS  =  A  : :  X0  =  X  : :  PRINT  "ANSWER" ;  CONS  : : 

PRINT  "X";  X0  ::  PRINT  ::  RETURN  ELSE  RETURN 
110  END 
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PROGRAM  TWO  (VARIABLE)  (See  Example  #10) 


X  =  ROW  ;  Y  =  COLUMN 

5  INPUT  "INIT  GRID?"  :  IG  : :  CONS  =  1.  E  +  100 

15  INPUT  "CROSSING"  :  LOOP  ::  INPUT  "N"  :  CUTTER  ::  INPUT 

"ROW  BEGINNING"  :  RB  : :  INPUT  "ROW  END"  :  RE  : :  INPUT 
"COLUMN  BEGINNING"  :  CB  :  :  INPUT  "COLUMN  END"  :  CE 
20  R1  =  RB  : :  GOSUB  80 

30  FOR  L123  =  1  TO  LOOP 

35  FOR  COL  =  CB  TO  CE  STEP  .5  A  CUTTER  ::  ANSWER  =  100* 

( (COL-R1  A  2)  2)  +  (6.4*  ( (COL  -  .5) A  2)  -  Rl  -  .6)  A  2  : : 
GOSUB  65  ::  NEXT  COL 

45  FOR  ROW  =  RB  TO  RE  STEP  . 5  A CUTTER  ::  ANSWER  =  100* 

((Cl  -  ROW  A  2)  A  2 )  +  (6 . 4*  (  (Cl  -  .5)/\2  -  ROW  -  .6)A2  :: 
GOSUB  70  : :  NEXT  ROW 
55  NEXT  L123 

60  GOTO  100 

65  IF  ANSWER  <  CONS  AND  Rl  >  COL  AND  COL  +  Rl  £  1  THEN 
CONS  =  ANS  ::  Cl  -  COL  ::  PRINT  "ANSWER"  ;  CONS  :: 

PRINT  "ROW"  ;  Rl  ::  PRINT  "COLUMN"  ;  Cl 

66  RETURN 

70  IF  ANSWER  <  CONS  AND  ROW  £  Cl  AND  Cl  +  ROW  £  1  THEN  CONS  = 
ANSWER  : :  Rl  =  ROW  : :  PRINT  "ANSWER"  :  CONS  : :  PRINT 
"ROW"  ;  Rl  ::  PRINT  "COLUMN"  ;  Cl 

71  RETURN 
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PROGRAM  TWO  (VARIABLE)  (Continued) 

75  IF  ANSWER  <  CONS  AND  ROW  >_  COL  AND  ROW  +  COL  £  1  THEN 
CONS  =  ANSWER  : : 

rL  =  row  :  :  Cl  =  COL  ::  PRINT  "ANSWER"  ;  CONS  ::  PRINT  "ROW"  ; 
Rl  ::  PRINT  "COLUMN"  :  Cl 

76  RETURN 

80  FOR  COL  -  CB  TO  STEP  (CE-CBJ/IG 
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PROGRAM  THREE  (VARIABLE)  (See  Example  #19) 


X  =  DEPTH;  Y  =  COLUMN  ;  2  =  ROW 

1  CONS  =  l.E  10 

10  CALL  CLEAR 

20  INPUT  "INITIAL  GRID"  :  IG 
30  INPUT  "ROW  START  &  END"  :  RB , RE 
40  INPUT  "COLUMN  START  &  END"  ;  CB,CE 

50  INPUT  "DEPTH  START  &  END"  :  DB ,  DE  ::  INPUT  "LOOP":L 
60  INPUT  "N"  :  N 

70  FOR  DEPTH  =  DB  TO  DE  STEP  (DE-DB) /IG 
80  FOR  COL  =  CB  TO  CE  STEP  (CE-CB)/IG 

90  FOR  ROW  =  RB  TO  RE  STEP  ( RE-RB) /IG  ::  A  =  ABS (DEPTH-1)  + 

ABS ( COL  -  1.5)  +  ABS (6*  ROW  -  1)  : :  GOSUB  500  ::  NEXT  ROW 

100  NEXT  COL  : :  NEXT  DEPTH 

105  PRINT  "OUT  OF  SUBPROGRAM" 

106  FOR  LI  =  1  TO  L 

110  FOR  DEPTH  =  DB  TO  DE  STEP  .5 -AN  ::  A  =  ABS  (DEPTH  -  1)  + 

ABS (COLO  -  1.5)  +  ABS (6*  ROWO  -  1)  GOSUB  600  ::  NEXT  DEPTH 

120  FOR  COL  =  CB  TO  CE  STEP  .5 -AN  ::  A  =  ABS  (DEPTHO  -  1)  + 

ABS (COL  -  1.5)  +  ABS (6*  ROWO  -  1)  ::  GOSUB  700  ::  NEXT  COL 
130  FOR  ROW  *  RB  TO  RE  STEP  .5/\  N  A  =  ABS  (DEPTHO  -  1)  + 

ABS (COLO  -  1.5)  +  ABS (6*  ROW  -  1)  ::  GOSUB  800  NEXT  ROW 

134  PRINT  LI 

135  NEXT  LI 


140  END 
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PROGRAM  THREE  (VARIABLE)  (Continued) 


500  IF  <  CONS  AND  DEPTH  +  COL  +2*  ROW  <_  3  THEN  DEPTHO  =  DEPTH 

: :  ROWO  =  ROW  : :  COLO  =  COL  : :  CONS  =  A  : :  GOSUB  1000  : : 

RETURN  ELSE  RETURN 

600  IF  A  <  CONS  AND  DEPTH  +  COLO  +  2*  ROWO  <  3  THEN  CONS  =  A  :: 

DEPTHO  =  DEPTH  ::  GOSUB  10 00  ::  RETURN  ELSE  RETURN 
700  IF  A  <  CONS  AND  DEPTHO  +  COL  +2*  ROWO  <_  3  THEN  COLO  =  COL 
::  CONS  =  A  GOSUB  1000  ::  RETURN  ELSE  RETURN 
800  IF  A  <  CONS  AND  DEPTHO  +  COLO  +2*  ROW  <_  3  THEN  CONS  =  A  :  : 

ROWO  =  ROW  : :  GOSUB  1000  : :  RETURN  ELSE  RETURN 

1000  PRINT  "ANSWER";  CONS:  "DEPTH" ? DEPTH : "COLUMN" ; COLO : "ROW" ; ROWO 
: :  PRINT  :  :  RETURN 
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COMPUTER  PRINT  OUT : RANDOM  SEARCH 


p 

00100 

C 

FIND  ABSOLUTE  MINIMUM  OF  NON-DIFFERENTABLE  BOUNED 

00200 

c 

FUNCTION  WITH  OR  WITHOUT  CONSTRAINT 

00300 

DIMENSION  A(3»50»50) t FMAX ( 50 ) >CMAX (50 » 3  >  fF  < 50 » 50 ) »CF<  50 

3)  f 

00400 

1 

MAX ( 50 1 50 )  » X  < 50 ) »Y(50) *  YMAX  <  50 » 50  > » CYMAX ( 50 » 50 »  3 ) fFMAXl 

50)  » 

00500 

2 

P<3>  fCMAXl (50 >3) f  CFAX < 50 1 3 ) >CL1 <50»50) » FCL ( 50 ) » CFCL ( 50 > 

00600 

3 

CMAX2 < 50 »  3 )  fFMAX2<50) » CL. <50*50*3)  »D(6) » 

00700 

4 

E  ( 3 ) 

00900 

N=20 

00900 

11=6 

01000 

12=3 

01100 

OPEN ( UNIT=20  »FILE= '  INI  ♦ DAT ' ) 

0 1  200 

OF'EN  <  UNI  T=21  » F ILE= '  RES « OUT  '  > 

01300 

I U. 1=20 

01400 

IU=5 

01500 

c 

READ  DOMAIN  PARAMETER 

01600 

READ  <IUI»#)(D<I)»I=1»6> 

*P 

01700 

01800 

5 

FORMAT ( 6F5 . 1 ) 

01900 

RAD=4 

02000 

N0UNT=0 

02100 

02200 

c 

INITALIZE  THE  ARRAY 

02300 

N=50 

02400 

DO  7  K=1 » I2»  1 

02500 

7 

A(K»1»1)=D<K) 

02600 

c 

CALCULATE  EPSILON  VALUES 

02700 

c 

02800 

DO  8  1  =  1 » 12  f 1 

02900 

8 

E ( I )  =  <  D ( 1+3 ) -D  < I > ) /50 • 0 

03000 

c 

03100 

c 

GENERATE  COORDINATES 

03200 

c 

♦P 

03300 

DO  40  J=1 fNf 1 

03400 

DO  30  1=1 fNf  1 

03500 

DO  9  K=1 » I2» 1 

03600 

9 

A<K»IfJ)=A<Kf1f1>*<~1>**<I*J*K>+I*E<K)*<-1)**I+J*E<K>*<- 

1) 

03700 

1 

**J+K*E(K)*<-1 )**N 

03800 

c 

03900 

c 

CHECK  IF  EXCEEDS  DOMAIN 

04000 

c 

04100 

DO  25  K= 1 > 1 2 » 1 

04200 

IF<A<KfIf J) .LE.D(K) )  A<Kf I f J)=D<K+3)-J*E(K> 

04300 

IF< A(Kf I f J) *GE .D<K+3) )  A(Kf I? J)=D<K)+J#E<K) 

04400 

25 

CONTINUE 

04500 

c 

04600 

c 

CALCULATE  THE  FUNCTIONAL  VALUES 

04700 

c 

p 


04900 

K1  =  I 

05000 

F  < I » J ) =CAL  <A»K1»N1> 

05.100 

N0UNT=K0UNT+1 

05200 

C 

WRITE! IUf 26) (A(K»I» J) »K=lf 12) ?  F  ( I » J ) 

05300 

26 

FORMAT (4X»3<E10»4f4X) t4X»E10.4/) 

05400 

30 

CONTINUE 

05500 

C 

05400 

C 

CALL  SUBROUTINES  TO  FIND  MIN  OF  THE  50  POINTS  JUST  EVALA 

TED . 

05700 

C 

05800 

K= J 

05900 

13=12 

06000 

CALL  FMAX I !  A  *  F » FMAX ,  CMAX  fKfNfI3) 

06100 

C 

06200 

c 

CALCULATION  TO  START  NEXT  SET  OF  POINTS 

06300 

c 

06400 

A  J= J 

*P 

06500 

AJ=AJ/50. 

06600 

DO  35  K=1»I2»1 

06700 

35 

A(K»1»1)=A(K»1»J) +A J*  < - 1 ) **  J 

06800 

40 

CONTINUE 

06900 

C 

07000 

C 

OUTPUT  MINIMUN  VALUE  8  COORDINATES 

07100 

C 

07200 

WRITE  !  IU »  45 ) 

07300 

45 

FORMAT <//»5X» 'VALUE  OF  THE  FUNCTION' »5Xf 'FIRST  COORDINAT 

E' f5Xf 

'SECOND 

07400 

1 COORD I NATE ' * 5X » ' THIRD  COORDINATE' ,/ ) 

07500 

DO  50  1=1 f N» 1 

07600 

50 

WRITE( IU*55)  FMAX !  I  >  f  ! CMAX ( I » K ) rK=l» 12) 

07700 

55 

F0RMAT!5XfE11.5f5XfE11.5f5XfEll,5f5XfE11.5) 

07800 

C 

07900 

C 

CALL  PLOTTING  ROUTINE 

08000 

C 

*P 

08100 

DO  200  J=1 f  4  f 1 

08200 

JJ=J 

08300 

CALL  MAX2<CMAX»FMAXf FMAX2fCMAX2»N» JJf 12) 

08400 

GO  TO  (60r62f6Af67)J 

08500 

60 

WRITE< IUf  61 ) J 

08600 

61 

FORMAT ( 10X » '  GROUP' f  13) 

08700 

WRITE! IUf 65)  FMAX2< J) t <CMAX2< JfK)fK=lfI2) 

08800 

65 

FORMAT  !5XfE12»4fE12.4fE12»4fE12*4f/) 

08900 

GO  TO  80 

09000 

62 

WRITE< IUf  61 ) J 

09100 

WRITE < IUf  65 )FMAX2! J ) »  <CMAX2<  J»K ) »K=1 » 12 ) 

09200 

GO  TO  80 

09300 

64 

WRITE ( IUf 61 )  J 

09400 

WRITE! IUf 65)  FMAX2 ( J  >  r  <  CMAX2 (JfK)?K=lfI2) 

09500 

GO  TO  80 

09600 

67 

WRITE! IUf 61 ) J 

* 
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P 


09700 

WRITE ( IU»  65 ) FMAX2  <  J ) » <  CMAX2 (J*K) »K=1>I2) 

09800 

80 

K=0 

09900 

DO  84  I = 1 » N  »  1 

10000 

IF <FMAX< I ) »EQ»9. 1E+10)GO  TO  84 

10100 

DIST=0 . 0 

10200 

DO  81  KI  =  1 » I2» 1 

10300 

P  <  K I  > =CMAX2  <  J » K I ) -CMAX ( I » K I ) 

10400 

81 

DIST=DIST+<P<KI )**2> 

10500 

Q=SGRT  <  DIST ) 

10600 

IF  <Q ♦ GE ♦ RAD )  GO  TO  84 

10700 

K*K+1 

10800 

FCL  <  K  >  =FMAX  < I ) 

10900 

DO  82  KI=1»I2»1 

11000 

82 

CFCL  <  K » K I ) =CMAX ( I , K I > 

11100 

FMAX (I)=9»1E+10 

11200 

*P 

DO  83  KI=1 » I2» 1 

11300 

83 

CMAX(I»KI)=9.1E+10 

11400 

84 

CONTINUE 

11500 

K=K+1 

11600 

FCL ( K ) =FMAX2 ( J ) 

11700 

DO  85  KI-1 » I2> 1 

11800 

85 

CFCL ( K  t K I ) =CMAX2  <  J , K I ) 

11900 

FMAX2< J)=9. 1E+10 

12000 

IF < K  .EO.  0.)  GO  TO  220 

12100 

GO  T0<101»lllfl21»131)J 

12200 

101 

K1=K 

12300 

DO  106  1=1 »K1  *  1 

12400 

CL1 ( J» I )=FCL< I ) 

12500 

FCL<I)*0. 

12600 

DO  105  KI=lfI2»l 

12700 

CL( J» I »KI )=CFCL< I >KI ) 

12800 

*P 

105 

CFCL( I »KI )*0. 

12900 

106 

CONTINUE 

13000 

GO  TO  200 

13100 

111 

K2=K 

13200 

DO  116  1*1 >K2» 1 

13300 

CL1 < J» I >=FCL< I ) 

13400 

FCL< I  >=0. 

13500 

13600 

DO  115  KI*1 » I2» 1 

13700 

CL(JrlfKI) *CFCL  < I » K I > 

13800 

115 

CFCL< I »KI )*0. 

13900 

116 

CONTINUE 

14000 

GO  TO  200 

14100 

121 

K3*K 

14200 

DO  126  1*1 »K3  f 1 

14300 

CLl(Jf I)*FCL(I) 

14400 

* 

FCL< I )*0* 

l 

l 
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p 

14500 

DO  125  KI  =  1 f 12  f 1 

14600 

CL ( J  f I f  K I ) =CFCL  ( I  f  K I ) 

14700 

125 

CFCL< I fKI >=0* 

14800 

126 

CONTINUE* 

14900 

GO  TO  200 

15000 

131 

K4=K 

15100 

DO  136  1=1 fK4f 1 

15200 

CL1 < J» I )=FCL< I ) 

15300 

FCL  < I ) =0 ♦ 

15400 

DO  135  KI=1fI2f1 

15500 

CL( J» I »KI )=CFCL< IfKI) 

15600 

135 

CFCL  < I »KI ) =0 . 

15700 

136 

CONTINUE 

15800 

200 

CONTINUE 

15900 

C 

16000 

C 

WRITE  THE  FUNCTION -VALUE  AND  THE  COORDINATES  OF  THE  UHOL 

E  CLUSTER 
*P 

16100 

c 

16200 

220 

IF(K1*EQ*0)  GOTO  225 

16300 

WRITE (IUf  222 )  K1 

16400 

222 

FORMAT  (IOXf'  CLUSTER  'fI3> 

16500 

DO  224  I=1fK1f1 

16600 

224 

WRITE (IUf*)  CL1(1fI)f (CL ( 1 f I fKI ) fKI=1 f 12 ) 

16700 

Ooer 

IF( K2 . EQ .0 )  GOTO  230 

16800 

WRITE < IUf226)K2 

16900 

226 

FORMAT (IOXf 'CLUSTER  'fI3> 

17000 

DO  228  1  =1 fK2  f 1 

17100 

228 

WRITE(IUf*)  CL1(2fI)f  <CL(2fIfKI) fKI  =  1fI2) 

17200 

230 

IF(K3.EG).0>  GO  TO  235 

17300 

WRITE( IUf232)K3 

17400 

232 

FORMAT ( IOXf 'CLUSTER  'fI3> 

17500 

DO  234  1  =  1 fK3  f 1 

17600 

*P 

234 

WRITE (IUf*)  CL1 (3 f I ) f (CL ( 3 f I fKI ) fKI=1 f 12 ) 

17700 

235 

IF(K4.EQ.O)  GO  TO  240 

17800 

WRITE (IUf 236)  K4 

17900 

236 

FORMAT (19Xf 'CLUSTER  'fI3) 

18000 

DO  238  1=1 fK4f l 

18100 

238 

WRITE (IUf*)  CL1 ( 4 f I ) f (CL (4 f I fKI ) fKI=1 f 12) 

18200 

240 

WRITE ( IUf  250) 

18300 
-  5Xf 

250 

FORMAT ( 5X f ' MAX  FUNCT  'fSXf'COORD  TE-1 ' f5Xf 'COORD  TE-2'f 

18400 

1 

'COORD  TE-3 ' f // ) 

18500 

KK=0 

.  18600 

DO  350  J=1 f  4  f 1 

18700 

GO  TO  (302f306f310f314)J 

18800 

302 

IF(KIfEQ.O)  GO  TO  350 

18900 

Jl»Kl 

19000 

I=J 

19100 

ICOUNT-O 

19200 

* 

303 

CALL  SMPLEX <  CL l f  CL  f  CYMAX  f  YMAX  f  D  f I f  J1 f 1 2  f ICOUNT ) 

X, 


p 

19300 

WRITE  < IU » 305 ) I COUNT 

19400 

305 

FORMAT <2X, 'GROUP  ITERATION' » 2X, 14 ) 

19500 

KK=KK+ICOUNT 

19600 

WRITE (IUf 304)  YMAX< J, 1 ) » <CYMAX( J» 1»KI)»KI=1»I2) 

19700 

304 

FORMAT (5X»E10*4»3<  6X  »E1 0 • 4 ) t //) 

19800 

GO  TO  350 

19900 

306 

IF (K2.EQ.0)  GO  TO  350 

20000 

J1=K2 

20100 

I=J 

20200 

GO  TO  303 

20300 

310 

IF(K3.EG.O)  GO  TO  350 

20400 

J1=K3 

20500 

I=J 

20600 

GO  TO  303 

20700 

314 

IF < K4 « EG « 0 )  GO  TO  350 

20800 

J1=N4 

*P 

20900 

I=J 

21000 

GO  TO  303 

21100 

350 

CONTINUE 

21200 

KNT-KK+KOUNT 

21300 

WRITE< IU»355)KNT 

21400 

355 

FORMAT (10X» 'TOTAL  ITERATION' ,2X»  15) 

21500 

CL0SE<UNIT=20*F.TLE='  INI  »DAT'  ) 

21600 

CLOSE ( UNIT=21 ,FILE= ' RES ♦ OUT ' ) 

21700 

STOP 

21800 

END 

21900 

SUBROUTINE  MAX2 ( CMAX3 » FMAX3 , FAX2 > CFAX2 ,  M ,  L  r  1 1 ) 

22000 

DIMENSION  FAX2(50) ,CMAX3<50,3) ,FMAX3(50> , CFAX2( 50 » 3 > ,TEM 

(3) 

22100 

11=3 

22200 

DO  10  1=1, M,1 

22300 

IF (FMAX3( I ) . GE.FMAX3< 1 ) )  GO  TO  10 

22400 

TEMP=FMAX3 ( 1 ) 

*P 

22500 

FMAX3( 1 )=FMAX3< I ) 

22600 

FMAX3 ( I ) =TEMP 

22700 

DO  5  K=1 , 3, 1 

22800 

TEM(K)=CMAX3( 1 »K) 

22900 

CMAX3  < 1 » K ) =CMAX3  < I , K ) 

23000 

CMAX3 ( I » K ) =TEM  <  K ) 

23100 

5 

CONTINUE 

23200 

10 

CONTINUE 

23300 

C 

23400 

C 

STORE  THE  CURRENT  LOWER  VALUE 

23500 

C 

23600 

FAX2(L)=FMAX3(1) 

23700 

FMAX3< 1)=9»1E+10 

23800 

DO  12  K=1 1 

23900 

CFAX2 <L»K> =CMAX3 ( 1 , K ) 

24000 

* 

12 

CMAX3<lfK)=9.1E+10 

24100 

RETURN 

24200 

END 

24300 

SUBROUTINE  FMAXI (B,H,FMAX1 ,CMAX1 ,L*N1 ,L3> 

24400 

DIMENSION  B ( 3 , 50 , 50 ) , H ( 50 » 50 ) » FMAXI (50) ,CMAX1 (50, 3) ,TEM( 

3) 

24500 

L3=3 

24600 

DO  10  1=1 ,N1 » 1 

24700 

IF(H(I,L)»GE.H(1,L)>  GO  TO  10 

24800 

TEMP=H( 1 ,L  > 

24900 

H(1,L)=H(I,L) 

25000 

H( I ,L )=TEMP 

25100 

DO  5  KK=1 , L3 , 1 

25200 

TEM(KK ) =B ( KK » 1 »L  > 

25300 

B(KK,1,L)=B(KK,I,L> 

25400 

5 

B(KK, I,L) =TEM(KK) 

25500 

10 

CONTINUE 

25600 

C 

*P 

25700 

C 

STORE  THE  CURRENT  LOWEST  VALUE  AND  COORDINATES 

25800 

c 

25900 

FMAX1 (L)=H(1»L) 

26000 

DO  12  KN=1 ,L3, 1 

26100 

12 

CMAX 1 ( L , KK  >  =B ( KK , 1 , L ) 

26200 

RETURN 

26300 

END 

26400 

C 

26500 

C 

LOCAL  OPTIMIZATION  USING  NELDER  AND  MEADS  SMPLEX  METHOD 

26600 

CC 

26700 

SUBROUT I NE  SMPLEX ( ACL 1 , CACL 1 , CZMAX , ZMAX , DD , L , M , I 4 , KOUT ) 

26800 

DIMENSION  ACL1 (50,50), CACL1 (50,50,3), CZMAX (50,50,3), ZMAX 

(50,50) 

, 

26900 

1  TEMPO ( 50 ) , DD ( 6 ) >  SUM2 (3),X2(4,4), SUM ( 3 ) 

27000 

DIMENSION  BST (4,4) »CBST (4,4,3) *CR( 4 , 3 ) »CP( 4 ,3 ) ,PIMG (4,3) 

27100 

1  FCR ( 4 ) »FPMG(4 ) , FCW( 4 ) »FEX ( 4 ) ,CW( 4 , 3 ) 

27200 

DIMENSION  XI (4) 

*P 

27300 

DIMENSION  EX (4, 3) 

27400 

RADD=0. 00004 

27500 

14=3 

27600 

ITER=0 ♦ 

27700 

N2=I4+1 

27800 

IU=5 

27900 

DO  3  1=1 »N2, 1 

28000 

BST(L , I)=+*91E+10 

28100 

DO  3  K=1 » 14, 1 

28200 

3 

CBST(L, I,K)=(+1)**K*0.91E+11 

28300 

WRITEdU,#)  BST  (L,1),BST(L,2),  BST  (L»3  >  ,BST  (L  ,4) 

28400 

C 

28500 

C 

FIND  THE  LOWEST  POINT  AND  THE  COORDINATES  BST(L,1) 

28600 

C 

28700 

WRITEdU, 5)  M 

28800 

5 

FORMA Td OX, 'THE  NUMBER  IS  ',14) 
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28900 

DO  10  1  =  1 f  M  f 1 

29000 

IF  < ACL1 (L  f  I  > .  EQ.0.91E+10) 

GOTO 

10 

29100 

IF < BST <L » 1 > ♦ LE . ACL 1(  L  *  I  > ) 

GO 

TO 

10 

29200 

TEMP=BST(L* 1) 

29300 

BST  ( L  f  1 ) =ACL 1( L  *  I ) 

29400 

ACLKLf  I  >=TEMP 

29500 

DO  8  K=1 *  14*1 

29600 

TEMPO (K)=CBST(L» 1* K) 

29700 

CBST<Lf lfK)=CACLl(Lf IfK> 

29800 

8 

CACL 1 ( L  f I f  K ) =TEMPO <  K  > 

29900 

10 

CONTINUE 

30000 

WRITE (IUf#)BST<Lfl> 

30100 

C 

30200 

C 

FIND  THE  SECOND  BEST 

30300 

C 

30400 

DO  20  1=1 f Mf 1 

*P 

30500 

IF ( ACL1 (LfI).EQ»BST(Lfl)  ) 

GO 

TO 

20 

30600 

IF (ACL1 <L f I )  .EQ,  0.91E+10) 

GOTO  20 

30700 

IF (BST  <Lf 2) .LE.ACL1 (Lf I ) >G0 

TO 

20 

30800 

TEMP=BST (Lf 2) 

30900 

BST ( L  f  2 ) =ACL 1 ( L  f I > 

31000 

ACLKLf  I)=TEMP 

31100 

DO  16  K=1 f 14  *  1 

31200 

TEMPO  <  K ) =CBST ( L  f  2  f  K ) 

31300 

CBST  <  L  f  2  f  K ) =CACL 1(L*I*K> 

31400 

16 

CACL1 (L  f I f K )=TEMPO<K ) 

31500 

31600 

20 

CONTINUE 

31700 

WRITE  < IU  f  # )  BST ( L  f  2 ) 

31800 

C 

31900 

C 

FIND  THE  THIRD  LOWEST  POINT 

32000 

C 

*P 

32100 

DO  26  1  =  1 f  M  f 1 

32200 

IF< ACLKLf  I)  .EG.BSKLf  1) ) 

GO 

TO 

26 

32300 

IF  (ACLKLf  I)  .ea.bst<lf2>> 

SIO 

to 

26 

32400 

IF ( ACL1 (Lfl) ♦EQ.0.91E+10) 

GOTO 

26 

32500 

IF  (BST  (Lf  3)  .LT.  ACLKLf  I) ) 

GO 

TO 

26 

32600 

TEMP=BST (L* 3) 

32700 

BST(Lf3)=ACLl(Lf I) 

32800 

ACLKLf  I  )=TEMP 

32900 

DO  25  K=1 f 14  f 1 

33000 

TEMPO (K)=CBST(L*3fK) 

33100 

CBST(Lf3fK)=CACLl(Lf IfK> 

33200 

25 

C ACL 1( L *  I  * K ) =TEMPO ( K > 

33300 

26 

CONTINUE 

33400 

WRITEdUf*)  BST (L * 3 ) 

33500 

c 

33600 

* 

C 

FIND  THE  FOURTH  LOWEST  POINT 

AT 

THIS  TIME 

p 


33700 

c 

33800 

DO  29  1=1 fMf 1 

33900 

IF ( ACL1 <LfI)*EQ»BST(Lf1>>  GO  TO  29 

34000 

IF ( ACL1  <  L  f  I )  *EQ . BST <!_ * 2 >  >  GO  TO  29 

34100 

IF<ACL1 (Lf  I ) . EQ.BST (L  *3) )  GO  TO  29 

34200 

IF < ACL! (LfI).GE.BST<Lf4))  GO  TO  29 

34300 

IF<ACL1<LfI>  .EQ.  0.91E+1 0)  GOTO  29 

34400 

TEMP=BST ( L  f  4 ) 

34500 

34600 

BST (Lf4)=ACL1 <Lf  I ) 

34700 

ACL1 <L  f  I ) =  TEMP 

34800 

DO  28  K=1 f  14  f  1 

34900 

TEMPO ( K ) =CBST  <  L  f  4  f  K ) 

35000 

CBST <  L  f  4  f  K ) =CACL KLfIfK) 

35100 

28 

CACL1 <L  » I fK)=TEMPO(K ) 

35200 

29 

CONTINUE 

*P 

35300 

WRITE ( IU» * )  BST < L f 4 ) 

35400 

C 

35500 

C  THE 

SIMPLEX  FORMED  BY  BST (L»  1  >  »BST (L»2)  rBST (Li»3) 

35600 

C  BST ( L * 4 >  WHICH 

35700 

35800 

C 

IS  THE  BIGGEST  ONE  SHOULD  BE  REMOVED 

35900 

WRITE <IUf 30)  BST(L»1) fBST(Lf2) »BST(Lf3) »BST(Lf4) 

36000 

30 

FORMAT (2X*F10.4»2X»F10. 4  f2X  fFIO  » 4»2X»F10»4) 

36100 

C 

36200 

C 

SEE  IF  THE  EXIT  CRITERIA  IS  SATISFIED 

36300 

C 

36400 

36500 

31 

DO  35  K=1 f 14 » 1 

36600 

SUM(K)=0. 

36700 

DO  35  1  =  1 fN2  f 1 

36800 

35 

SUM<K)=SUM(K)+CBST(LfIfK) 

*P 

36900 

DO  36  k  =  l f 14  f 1 

37000 

36 

XI  (K)=SUM(K')/4  ♦ 

37100 

DIF=0 . 

37200 

DO  40  1=1 fN2f 1 

37300 

DO  39  K2=l f 14 f 1 

37400 

X2 (If  K2 ) =CBST  <L  f I f  K2 ) -X 1 ( K2 ) 

37500 

39 

dif=dif+x2< i fK2)*#2 

37600 

40 

CONTINUE 

37700 

DIFB=SQRT (DIF) 

37800 

ITER=ITER+1 

37900 

IF<DIFB.LE.RADD)  GO  TO  500 

38000 

URITE<IUf43)  DIFB 

38100 

43 

FORMAT < 4X f ' DIFFERENCE*  'fE12.5) 

38200 

C 

REMOVE  THE  HIGHEST  FROM  THE  SIMPLEX 

38300 

C 

THE  HIGHEST  POINT  IS  THE  BST(Lf4) 

38400 

C 

THE  MEDIAN  OF  THE  POINTS  BST<L f 1 ) f BST <L f 2) f BST(L f 3 ) 

* 
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38500 

DO  48  K=1 » 14 » 1 

38600 

SUM2<K)=0. 

38700 

N3=N2-i 

38800 

DO  46  1=1 fN3f 1 

38900 

46 

SUM2 <  K ) =SUM2  <  K ) +CBST <  L  f  I  f  K  > 

39000 

48 

CP <  L  f  K  > =SUM2  <  K ) /3 . 

39100 

C 

39200 

c 

FIND  THE  IMAGE  OF  BST(L*4)  THOUGH  CP 

39300 

c 

39400 

DO  50  K=1 » 14 » 1 

39500 

50 

PIMG<LfK)=2.#CP<LfK)-CBST (L*4»K> 

39600 

c 

39700 

c 

CHECK  IF  EXCEEDS  THE  DOMAIN  OF  THE  FUNCTION 

39800 

c 

39900 

DO  52  K=1 f  14  f  1 

40000 

I F ( P IMG ( L  f  K ) . LT . DD ( K ) )  P IMG  <  L » K ) =DD  < K ) 

*P 

40100 

52 

IF<PIMG<LfK) *GT . DD ( K+3 ) )  PIMG<LfK)=DD<K+3> 

40200 

c 

40300 

c 

EVALUATE  THE  FUCTION  AT  THESE  POINTS 

40400 

c 

40500 

L3=L 

40600 

M1  =  I4 

40700 

FPMG(L)=XFCT  (PIMGfMI fL3) 

40800 

K0UT=K0UT+1 

40900 

write(iu»54)  fproa<l) 

41000 

54 

FORMAT < 4X » /FPMG=/ fE12.4) 

41100 

IF (FPMG(L )  .  LT , BST  <L f 1 ) )  GO  TO  200 

41200 

IF (FPMG<L) .LT .BST <Lf2)  )  GO  TO  100 

41300 

IF (FPMG (L ) . GT . BST <L  f  4 ) )  GO  TO  60 

41400 

DO  56  K=1 fI4f1 

41500 

56 

CR(LfK)=<3*CP<LfK)-CBST<Lf4fK> >/2. 

41600 

c 

*P 

41700 

c 

CHECK  IF  EXCEEDS  DOMAIN 

41800 

c 

41900 

DO  58  K=1 fK4f  1 

42000 

IF<CR<LfK> .LT.DD(K) )  CR(LfK)=DD(K) 

42100 

58 

IF<CR  <L  fK )  .GT.  DD(K+3) )  CR(L»K>=DD(K+3) 

42200 

42300 

c 

42400 

c 

EVALURATE  THE  FUNCTION 

42500 

c 

42600 

M1  =  I4 

42700 

L3=L 

42800 

FCR<L>=XFCT (CRfMI fL3) 

42900 

K0UT=K0UT+1 

43000 

BST  <Lf4)=FCR<L) 

43100 

DO  59  K=1 f  14  f  1 

43200 

59 

CBST<Lf4fK)-CR<LfK> 

* 


.  I  I  A •  in  Him  n  M  i  m  j 
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43300 

GO  TO  400 

43400 

60 

DO  65  K=1 f 14  f 1 

43500 

65 

CW<LfK>=<CP<LfK)+CBST<Lf4fK> >/2. 

43600 

C 

43700 

C 

SEE  IF  EXCEEDS  DOMAIN 

43800 

C 

43900 

DO  70  K=1 f I4f 1 

44000 

IF<CW<LfK>  .LT.  DD<K>>  CU<LrK>=DIKK> 

44100 

44200 

70 

IF(CW(LfK) .GT.  DD<K+3 ) )  CW(LfK)=DD<K+3> 

44300 

M1  =  I4 

44400 

L3=L 

44500 

FCW(L)=XFCT (CWfM1fL3> 

44600 

K0UT=K0UT+1 

44700 

BST<Lf4)=FCW<L> 

44800 

*P 

DO  75  K=1 f  14  f  1 

44900 

75 

CBST <  L  f  4  f  K ) =CW  <  L » K ) 

45000 

GO  TO  400 

45100 

100 

BST  <  L » 4 ) =FPMG  <  L ) 

45200 

DO  110  K-l f  14  f  1 

45300 

110 

CBST (Lf4fK)=PIMG(LfK) 

45400 

45500 

GO  TO  400 

45600 

200 

DO  220  K=1 » 14  f 1 

45700 

220 

EX  <  L  f  K ) =3 « 0#CP  <  L  f  K ) -2 . 0#CBST  <  L  f  4  f  K ) 

45800 

C 

45900 

C 

SEE  IF  EXCEEDS  DOMAIN 

46000 

C 

46100 

DO  250  K5=l f  14  f  1 

46200 

IF(EX(LfK5) .LT.  DD<K5) )EX<LfK5>=DD(K5+3) 

46300 

46400 

*P 

250 

IF(EX(LfK5) .GT.DD<K5+3) )  EX<LfK5)=DD(K5+3) 

46500 

M1  =  I4 

46600 

L3=L 

46700 

FEX<L)=XFCT(EXfM1»L3) 

46800 

K0UT=K0UT+1 

46900 

IF<FEX<L> .LT.BST<Lf1> >  GO  TO  310 

47000 

WRITE<IUf255>  FEX<L> 

47100 

255 

FORMAT < 1 OX f 'FEX= ' fEIO  .  4 ) 

47200 

GO  TO  100 

47300 

310 

BST  <  L  f  4 ) *FEX  <  L ) 

47400 

DO  320  K=1fI4f1 

47500 

320 

CBST  <  L  f  4  f  K ) =EX  <  L  f  K ) 

47600 

400 

DO  410  1=1 »N2f 1 

47700 

IF < BST <Lf 1 ) .LT .BST (Lf I ) )  GO  TO  410 

47800 

TEMP=BST<Lf1> 

47900 

BST(Lf1)*BST<LfI> 

48000 

* 

BST<LfI>=TEMP 
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48100 

DO  405  K=1 f 14  f 1 

48200 

TEMPO (K ) =CBST  <  L  f 1 fK ) 

48300 

CBST  <  L  f  1  f  K ) =CBST ( L  f  I  f  K ) 

48400 

405 

CBST  <  L » I » K ) =TEMPO ( K ) 

48500 

410 

CONTINUE 

48600 

DO  450  1=2  f  N2f 1 

48700 

IF< BST (L » 2) »LT .BST (L » I ) )  GO  TO  450 

48800 

TEMP=BST  <  L  *  2 ) 

48900 

BST <Lf2)=BST <Lf  I ) 

49000 

BST  (L  f  I ) =TEMP 

49100 

DO  420  K=1 f  14  f  1 

49200 

TEMPO ( K ) =CBST <Lf2fK) 

49300 

CBST ( L.  f  2  f  K  ) =CBST <  L  f  I  f  K  > 

49400 

420 

CBST ( L » I f  K ) =TEMPO ( K ) 

49500 

450 

CONTINUE 

49600 

IF  (BST  (L.  r  3 )  .LT  .  BST  <L  f  4  )  )  GO  TO  460 

*P 

49700 

TEMP=BST(Lf3> 

49800 

BST  <L  f 3>=BST (L.  f  4) 

49900 

BST  <Lf4)=TEMP 

50000 

DO  455  K=1 » 14  f 1 

50100 

TEMPO  <  K ) =CBST  <  L  f  3  f  K ) 

50200 

CBST <Lf3fK)=CBST (Lf4fK) 

50300 

455 

CBST  <  L  f  4  f  K ) =TEMPO  <  K ) 

50400 

460 

IF< ITER.LT ,10. )GOTO  31 

50500 

ITER=0. 

50600 

BST <L » 4 )  =  <BST <L » 1 )+BST <IL  f  2>+BST (L  f3>+BST (Lf4>>/4. 

50700 

GO  TO  400 

50800 

500 

ZMAX  <  L » 1 ) =BST  <  L  f  1 ) 

50900 

DO  505  K=1 f  14  f  1 

51000 

505 

CZMAX ( L » 1 » K ) =CBST  <  L » 1 »  K ) 

51100 

WRITE <IU .510)  ZMAX(LfI) » (CZMAX(L f 1 fK) fK=1 f 14 > 

51200 

510 

FORMAT ( 10XfE10.4f3(E10.4f4X) f/> 

*P 

51300 

RETURN 

51400 

END 

51500 

FUNCTION  XFCT(CfIIfIP) 

51600 

DIMENSION  C(4f3) 

51700 

11=3 

51800 

XFCT=9 .0-8 ,0#C <  IP f  1 )  -6 . 0*C <  IP  f2 )-4 .0#C ( I 

51900 

1 

Pf3)+2.0#C<IPf1)**2+2.0*C(IPf2>**2+C(IPf3>*#2 

52000 

2 

fC( IPf 1 )#C< IP f 2 ) *2 . 0+2 . 0#C < IPf1)*C(IPf3) 

52100 

RETURN 

52200 

END 

52300 

FUNCTION  CAL(C1fL3fN3) 

52400 

DIMENSION  Cl <3f50f50) 

52500 

CAL=9.0-8.0*C1(1fL3fN3)-6.0*C1<2fL3fN3)-4.0*C1(3fL3fN3)+ 

52600 

1 

2.0*C1<1fL3fN3>**2+2.0*C1(2fL3fN3>**2+C1<3fL3fN3>**2 

52700 

2 

+2.0*C1(1fL3fN3)#C1<2fL3fN3>+2.0*C1<1fL3fN3>*C1<3fL3fN3> 

52800 

RETURN 

* 


